PL-TR-97-2012 


APPLICATION  OF  PARALLEL  COMPUTING 
ALGORITHMS  TO  THREE-DIMENSIONAL  MODELING 
OF  SEISMIC  WAVES  GENERATED  BY  COMPLEX 
SOURCES  IN  INHOMOGENEOUS  MEDIA 


Charles  G.  Doll,  Jr. 
Joseph  R.  Matarese 
Bertram  Nolte 


Earth  Resources  Laboratory 
Department  of  Earth,  Atmospheric,  and 
Planetary  Sciences 

Massachusetts  Institute  of  Technology 
42  Carleton  Street 
Cambridge,  MA  02142 


15  January  1997 


Final  Technical  Report 
14  July  1995-30  November  1996 


19970606  131 


Approved  for  Public  Release;  Distribution  Unlimited 


PHILLIPS  LABORATORY 
Directorate  of  Geophysics 
AIR  FORCE  MATERIEL  COMMAND 
HANSCOM  AFB,  MA  01731-3010 


lync  QUALTpy  ir;sPECTED  i 


The  views  and  conclusions  contained  in  this  document  are  those  of  the  authors  and 
should  not  be  interpreted  as  representing  the  official  policies,  either  express  or  implied,  of 
the  Air  Force  or  U.S.  Government. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


'jAMES  F.  LEWKOWI^ 
Alternate  Contract  Manager 
Earth  Sciences  Division 


/ 

fAMES  F.  LEWKOWICj/ 
Director 

Earth  Sciences  Division 


This  report  has  been  reviewed  by  the  ESD  Public  Affairs  Office  (PA)  and  is  releasable  to 
the  National  Technical  Information  Service  (NTIS). 


Qualified  requestors  may  obtain  copies  from  the  Defense  Technical  Information  Center. 
All  others  should  apply  to  the  National  Technical  Information  Service. 


If  your  address  has  changed,  or  you  wish  to  be  removed  from  the  mailing  list,  or  if  the 
addressee  is  no  longer  employed  by  your  organization,  please  notify  PL/IM,  29  Randolph 
Road,  Hanscom  AFB,  MA  01731-3010.  This  will  assist  us  in  maintaining  a  current 
mailing  list. 


Do  not  return  copies  of  this  report  unless  contractual  obligations  or  notices  on  a  specific 
document  requires  that  it  be  returned. 


REPORT  DOCUMENTATION  PAGE 


Form  Approved 
0MB  No.  0704-0188 


Publfc  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruaions,  searching  existing  data  sources 
gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  colleaion  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  asoect  of  this 
collection  of  information,  including  ^^^gfestions  for  reducmg  this  bidden  to  Washington  Headquarters  Services.  Oireaorate  for  Information  Operations  and  Reacts,  1215  Jefferson 
Oavis  Highway,  Suite  1204,  Arlington,  VA  22202-4302.  and  to  the  Office  of  Management  and  Budget.  Paperwork  Reduction  Projea  (0704-0188),  Washington.  DC  20503. 


2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

15  JAN  1997  Final  Technical  7/14/95-11/30/96 


4.  TITLE  AND  SUBTITLE 


1.  AGENCY  USE  ONLY  (Leave  blank) 


^plication  of  Parallel  Computing  Algorithms  to. Three- 
Dimensional  Modeling  of  Seismic  Waves  Generated  by 
Complex  Sources  in  Inhomogeneous  Media 


6.  AUTHOR(S) 

Charles  G.  Doll,  Jr.  . 
Joseph  R.  Matarese 
.Bertram  Nolte 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Earth  Resources  Laboratory 

Dept,  of  Earth,  Atmospheric,  and  Planetary  Sciences 
Massachusetts  Institute  of  Technology 
42  Carleton  Street 
Cambridae,  MA  02142 


9.  SPONSORING /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

Phillips  Laboratory 
29  Randolph  Road 
Hanscom  AFB,  MA  01731-3010 

Contract  Manager:  Delaine  Reiter/GPE 


5.  FUNDING  NUM3ERS 


F19628-95-K-0008 


PE  62101F 
PR  7600 


TA  GM 


WU  AA 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


63614 


10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 


PL-TR-^7-2012 


12a.  DISTRIBUTION /AVAILABILITY  STATEMENT 


12b.  DISTRIBUTION  CODE 


Approved  for  Public  Release;  Distribution  Unlimited 


13.  ABSTRACT  (Maximum  200  words)  ^  - 

We  describe  a  hardware  upgrade  that  was  performed  on  the  nCUBE  2 
massively  parallel  computer  at  MIT*s  Earth  Resources  Laboratory. 

The  upgrade  involved  the  installation  of  additional  memory,  additional 
disk  space,  and  an  additional  small  nCUBE  2  computer.  The  reason  for 
performing  the  upgrade  was  to  provide  ERL  with  a  computational  environment 
suitable  for  conducting  large-scale  2-D  and  3-D  wave  propagation 
simulations.  We  also  describe  two  recent  projects  that  were  performed 
on  the  upgraded  nCUBE  2.  One  project  was  a  3-D  simulation  of  earthquake¬ 
generated  ground  motion  in  the  Boston  Basin.  The  other  project  was  the 
development  of  a  parallel  irregular-grid  modeling  code  which  is  well-suited 
for  2-D  modeling  of  regional  wave  propagation  in  the  presence  of  surface 
topography . 


14.  SUBJECT  TERMS  ^  ~ 

Parallel  computing;  3-D  seismic  wave  propagation; 
earthquake  ground  motion 


15.  NUMBER  Or  PAGES 

20 _ 

16.  PRICE  CODE 


17.  SECURITY  CLASSIFICATION  118.  SECURITY  CLASSIFICATION  Il9.  SECURITY  CLASSIFICATION  20.  LIMITATION  OF  ABSTRACT  • 
OF  REPORT  I  OF  THIS  PAGE  I  OF  ABSTRACT 

Unclassified  I  Unclassified  |  Unclassified  SAR 


NSN  7540-01*280-5500 


Standard  298  (Rev  2-89) 

Pr fbCr-DC'd  0'/ 


TABLE  OF  CONTENTS 


List  of  Contributing  Scientists  . 

List  of  Previous  and  Related  Contracts . 

Bibliography  of  Publications  Totally  or  Partially  Supported  by  the  Contract 

Introduction  . 

Boston  Basin  Ground  Motion  Simulations . 

Irregular-Grid  Modeling  . 

References  . 


List  of  Contributing  Scientists 

Charles  G.  Doll,  Jr.,  Research  Specialist,  Massachusetts  Institute  of  Technology. 
Joseph  R.  Matarese,  Visiting  Scientist,  Massachusetts  Institute  of  Technology. 
Bertram  Nolte,  Postdoctoral  Associate,  Massachusetts  Institute  of  Technology. 


List  of  Previous  and  Related  Contracts 

DARPA/AFGL  Contract  F19628-89-K-0020,  “Regional  Seismograms:  Attenuation  and  Scat¬ 
tering”,  July  1989  to  June  1991. 

DARPA/AFPL  Contract  F19628-90-K-0057,  “Research  in  Regional  Seismology:  The  Effect 
of  Anisotropy”,  August  1990  to  July  1992. 

DARPA/AFPL  Contract  F29601-91-K-DB15,  “Research  on  Monitoring  at  Regional  Dis¬ 
tances”  ,  September  1991  to  September  1993. 

AFOSR  Contract  F49620-92-J-0413,  “Basic  Research  in  Nuclear  Test  Monitoring”,  July  1992 
to  July  1994. 

AFOSR  Contract  F49620-93- 1-0424  “Seismic  Wave  Radiation,  Propagation  and  Event  Lo¬ 
cation  in  Laterally  Heterogeneous  Media”,  July  1993  to  December  1994. 

AFOSR  Grant  F49620-94-1-0282,  “Effect  of  3D  Heterogeneities  and  Topography  on  Seismic 
Wave  Propagation  and  the  Use  of  Empirical  Green’s  Functions  for  Source  Characteriza¬ 
tion  and  Discrimination”,  May  1994  to  May  1996. 

AFOSR  Grant  F49620-94-1-0273,  “Characterization  of  Seismic  Sources  Using  Empirical 
Green’s  Functions”,  July  1994  to  June  1997. 

AFOSR  Grant  F19628-95-C-0091,  “Research  on  Seismic  Monitoring  at  Regional  Distances,” 
June  1995  to  June  1996. 

Bibliography  of  Publications  Totally  or  Partially  Sponsored  by 

the  Contract 

Olsen,  K.B.,  R.J.  Archuleta,  and  J.R.  Matarese,  1995,  Three-dimensional  simulation  of  a 
magnitude  7.75  earthquake  on  the  San  Andreas  Fault,  Science,  270,  1628-1632. 


1.  INTRODUCTION 


The  512-processor  nCUBE  2  massively  parallel  computer,  originally  purchased  in  1990,  was 
upgraded  to  provide  a  capable  computing  environment  for  conducting  large-scale  2-D  and 
3-D  wave  propagation  simulations  at  the  M.I.T.  Earth  Resources  Laboratory  (ERL).  This 
upgrade  involved  three  major  components:  (1)  additional  memory  to  allow  simulations  to 
be  conducted  on  larger  earth  models;  (2)  additional  disk  space  to  allow  detailed  results  of 
these  simulations  to  be  recorded  for  further  analysis;  and  (3)  an  additional,  small  nCUBE  2 
computer  to  be  used  ultimately  as  an  Oracle  database  platform  and  to  enable  code  develop¬ 
ment  work  to  continue  while  large-scale  simulations  were  occupying  the  big  machine.  The 
memory  upgrade  was  performed  in  late  December,  1995,  and  the  disk  upgrade  and  secondary 
nCUBE  2  installation  were  performed  in  the  early  spring  of  1996. 

Quantitatively,  the  memory  upgrade  increased  the  nCUBE  2’s  aggregate  memory  from 
2.75  GB  (gigabytes)  to  5.375  GB,  while  the  disk  upgrade  added  a  100  GB  high  performance 
parallel  filesystem.  Before  the  memory  upgrade,  the  nCUBE  2  had  448  processors  with  4 
MB  (megabytes)  of  RAM  each  and  ,64  processors  with  16  MB  of  RAM  each.  Given  this 
configuration  and  the  nature  of  the  parallel  finite  difference  modeling  algorithm  used  in  the 
wave  simulations,  the  largest  models  that  could  be  stored  in  memory,  decomposed  among  all 
512  processors,  comprised  500  by  500  by  180  grid  points.  The  upgrade  has  quadrupled  the 
memory  on  roughly  one-half  the  machine,  i.e.,  the  memory  on  224  processors  was  increased 
from  4  MB  to  16  MB.  Therefore,  the  maximum  model  size  has  approximately  doubled,  for 


1 


example,  a  1000  by  500  by  180  grid  can  now  be  stored. 

In  the  course  of  conducting  finite  difference  simulations,  snapshots  of  the  surface  particle 
velocity  field  are  typically  recorded  at  every  tenth-time  step  to  provide  a  record  of  the  surface 
ground  motion.  Given  the  duration  of  the  simulations,  the  files  containing  the  snapshot  data 
grow  to  several  hundred  megabytes  of  data  for  each  of  the  three  Cartesian  components  of  the 
particle  velocity.  The  size  of  the  disk  array  now  makes  it  possible  to  archive  the  results  of 
many  simulations  for  comparative  analysis.  Also,  snapshots  of  particle  motion  for  the  entire 
model  volume  can  be  saved,  although  at  100  MB  per  snapshot,  not  at  the  same  rate  as  with 
the  surface  motion. 

The  following  sections  describe  examples  of  wave  propagation  simulation  performed  on 
the  upgraded  nCUBE  2  computer. 

2.  BOSTON  BASIN  GROUND  MOTION  SIMULATIONS 

In  a  previous  study  (Olsen  et  al.,  1995),  a  parallel  version  of  a  3-D  elastic  finite-difference 
code  was  run  on  ERL’s  nCUBE  2  512-processor  computer  to  simulate  earthquake-generated 
ground  motions  in  the  Los  Angeles  Basin.  The  purpose  was  to  study  site  amplification 
and  scattering  effects  due  to  the  3-D  structure  of  the  basin.  Here  we  report  preliminary 
results  of  a  similar  study  of  the  Boston  Basin.  Using  the  same  finite  difference  code  on  our 
nCUBE  2,  we  simulated  the  3-D  ground  motions  within  the  Boston  Basin  generated  by  a 
hypothetical  earthquake  in  the  basin.  The  results  of  some  initial  simulations  for  the  Boston 
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Basin  suggested  that  the  3-D  structure  of  the  sedimentary  basin  has  a  significant  effect  on 
the  complexity  and  amplification  of  surface  ground  motions.  To  better  understand  the  basin 
effects,  we  have  done  subsequent  simulations,  described  below,  that  compare  the  ground 
motions  resulting  from  the  3-D  basin  model  with  those  from  a  simple,  1-D  model  (layer  over 
a  half-space)  that  does  not  contain  a  3-D  sedimentary  basin. 

For  each  earth  structure  model,  seismic  waves  were  simulated  for  a  hypothetical  Mw  5.8 
earthquake.  The  earthquake  occurs  along  a  vertical,  strike-slip  fault  striking  east- west  with 
length  6.0  km  and  height  3.8  km.  The  fault  centroid  (point  of  initiation  of  the  rupture)  is 
placed  at  about  a  5  km  depth,  and  the  fault  was  allowed  to  displace  1  meter  to  the  west. 
For  the  simulation  done  with  the  3-D  basin  model,  the  top  of  the  fault  lies  400  meters  closer 
to  the  surface  than  for  the  1-D  model. 

Each  earth  structure  model  (3-D  and  1-D)  was  gridded  with  a  cell  size  of  100  meters, 
which  permits  adequate  sampling  of  wavelengths  down  to  ~600  meters.  Given  the  P  propa¬ 
gation  velocity  of  ~5000  meters/sec  in  the  sedimentary  rock  of  the  basin,  the  peak  frequency 
of  the  simulated  ground  motion  is  8  Hz,  with  a  center  frequency  of  4  Hz.  The  size  of  the  basin 
model  grid  was  401  by  400  by  72  nodes  in  the  east-west,  north-south  and  depth  directions, 
respectively. 

The  results  of  the  ground  motion  simulations  are  shown  in  Figure  1  (3-D  basin  model) 
and  Figure  2  (1-D  two-layer  model).  From  the  top  frames  of  the  figures  we  see  that  the  3-D 
basin  model  produces  peak  ground  velocities  that  are  strongly  variable  with  position  on  the 
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BOSTON  BASIN 
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Figure  1:  Snapshot  of  peak  ground  velocity  at  the  surface  of  the  Boston  Basin  resulting 
from  an  Mw  =  5.8  earthquake  (top).  Ground  motion  at  surface  25  seconds  after  start  of 
rupture  (bottom).  The  source  of  the  earthquake  is  a  vertical  strike-slip  fault  displaced 
1  m  to  the  west.  The  fault  has  a  length  of  6  km,  a  height  of  3.8  km,  a  centroid  about 
5  km  deep,  and  is  located  at  the  southeast  edge  of  the  basin  (note  its  lobed  radiation 
pattern  at  the  southeast  corner  of  the  basin).  The  earth  model  consists  of  a  3-D  basin 
of  sedimentary  rock  and  irregular  geometry  surrounded  by  hard  rock  (granite). 
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LAYER  OVER  HALFSPACE 
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Figure  2:  Peak  ground  velocity  at  the  surface  using  a  layer-over-a-half-space  earth  structure 
model  (top).  Snapshot  of  ground  motion  taken  25  seconds  after  rupture  initiation  (bot¬ 
tom).  Source  parameters  of  the  hypothetical  earthquake  and  overall  dimensions  of  the 
earth  model  are  the  same  as  for  the  basin  example  in  Figure  1.  The  layer-over-a-half-space 
model  consists  of  a  2  km  thick  sedimentary  rock  layer  over  a  hard  rock  (granite)  half¬ 
space.  Note  the  spatial  uniformity  of  the  peak  ground  velocity  (top)  and  the  simplicity 
of  the  ground  motion  (bottom). 


earth’s  surface.  In  contrast,  peak  velocities  for  the  1-D  model  (Figure  2)  are  rather  uniform. 
Also,  on  average,  the  1-D  model  yields  lower  peak  ground  velocities.  For  the  3-D  model 
(Figure  1),  the  high  peak  ground  velocities  forming  the  lobes  of  the  radiation  pattern  in  the 
immediate  vicinity  of  the  fault  (southeast  corner  of  the  basin)  is  explained  by  the  fact  that 
the  top  of  the  fault  is  shallower  than  in  the  1-D  model. 

The  bottom  frame  of  Figures  1  and  2  compare  snapshots  of  the  surface  ground  motion 
for  the  two  models.  Each  snapshot  corresponds  to  the  wavefield  25  seconds  after  rupture 
initiation.  We  see  that  the  complex  pattern  of  ground  motion  associated  with  the  basin 
model  (bottom  frame  of  Figure  1)  stands  in  sharp  contrast  to  the  simple  wavefield  in  the 
layer-over-a-half-space  model  (Figure  2).  The  wavefield  in  the  1-D  model  is  attributed  to 
the  primary  wavefront  radiating  outward  from  the  source,  while  the  wavefield  in  the  3-D 
basin  model  is  a  superposition  of  the  primary  wavefront  with  waves  scattered  by  the  basin 
structure. 

Suggested  future  work  includes  expanding  the  volume  of  the  basin  model  to  ensure  against 
possible  contamination  from  artificial  ground  motion  reflections  generated  at  the  absorbing 
boundaries.  Further,  a  major  modification  of  the  finite  difference  code  is  the  addition  of 
anelastic  attenuation  operators.  The  interaction  of  anelastic  attenuation  with  the  effects 
of  basin  geometry  (focusing  and  scattering  of  seismic  waves)  will  yield  a  realistic  estimate 
of  the  absolute  levels  of  ground  motion  and  refine  the  extent  of  ground  motion  variability 
associated  with  a  sedimentary  basin  structure. 
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3.  IRREGULAR-GRID  MODELING 


Another  research  project  was  the  development  of  an  irregnlar-grid  method,  based  on  a  finite 
volume  approximation,  for  the  modeling  of  wave  propagation  in  crustal  models  containing 
irregular  interfaces  and  surface  topography.  This  method  has  the  advantage  over  finite  dif¬ 
ferences  in  that  interfaces  in  a  model  (including  an  irregular  free  surface)  can  be  represented 
more  accurately,  and  that  the  grid  spacing  can  vary  throughout  the  model  and  can  thus  be 
scaled  to  the  wavespeeds. 

We  have  implemented  2-D  parallel  versions  of  the  algorithm  for  both  SH  and  P-SV 
waves  on  the  nCUBE  2  parallel  computer  of  the  Earth  Resources  Laboratory,  using  the 
Message  Passing  Interface  (MPI).  Using  MPI  guarantees  easy  portability  of  the  codes. 
An  efficient  parallel  implementation  of  the  irregular-grid  method  was  considerably  more 
challenging  than  a  parallel  implementation  of  a  finite  difference  algorithm,  as  an  efficient 
grid  decomposition  (in  terms  of  both  load  balancing  and  computation-to-communication 
ratio)  is  far  more  complicated  for  an  unstructured  grid. 

We  used  the  irregular-grid  code  to  study  the  effect  of  surface  topography  on  regional 
wave  propagation.  Figure  3a  shows  a  2-D  crustal  model  with  surface  topography.  Figure 
3b  shows  the  topography  in  more  detail.  The  topography  was  obtained  with  the  online 
profile  maker  from  Cornell  University’s  Middle  East  and  North  Africa  Project  database 
{http://atlas.geo.comell.edu)  for  a  profile  in  northern  Iran. 

We  computed  synthetic  seismograms  for  this  model  with  and  without  including  the  sur- 
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face  topography.  The  results  are  shown  in  Figure  4.  There  are  some  distinct  differences  in 
the  waveforms  obtained  for  the  two  cases.  In  this  example  the  most  prominent  differences 
are  on  the  vertical  component  at  larger  arrival  times  (greater  than  70  s).  More  arrivals  with 
higher  amplitudes  are  present  in  the  case  of  surface  topography.  These  arrivals  are  most 
likely  due  to  surface  waves  that  arise  from  scattering  at  the  rough  surface. 

Our  comparison  demonstrates  that  a  pronounced  surface  topography,  such  as  is  present  in 
many  regions  (e.g.,  the  Middle  East),  will  have  a  distinct  effect  on  observed  waveforms,  and 
should  thus  be  included  in  the  modeling.  Our  newly  developed  method  provides  a  valuable 
tool  for  this  task. 

Suggested  future  work  in  the  area  of  irregular-grid  modeling  includes  the  extension  of 
the  method  to  3-D  and  its  coupling  to  regular-grid  finite  differences,  which  will  enable  us 
to  use  irregular  grids  only  in  parts  of  the  model,  particular  in  the  near-surface  region.  Such 
a  hybrid  approach  would  take  advantage  of  the  strong  points  of  either  method  and  thus 


maximize  the  overall  efficiency. 
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Figure  4:  Synthetic  waveforms  on  horizontal  and  vertical  components  for  the  model  in 
Figure  3  and  for  an  equivalent  model  with  a  flat  surface. 
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